pacman::p_load(dplyr,tidyr, rcompanion, readxl,ggplot2, ggstatsplot,brms, bayesplot, rstan, tidybayes,ggdist)
df<-read_excel("Food_Habits_109_sex_answers.xlsx")%>%
dplyr::select(`Group (as in the html file: 1, CoV_smell_taste_disturbances; 2, CoV_no_smell_taste_disturbances; 3, healthCtr)`,Number, Sex,
`How much salt do you use in your  meals?`,
`How much bitter additives do you use in your  meals?`,
`How much sugar/sweeteners do you use in your meals?`,
`How much sour additives (e.g. citric acid, vinegar) do you use in your  meals?`)%>%
rename("Group"="Group (as in the html file: 1, CoV_smell_taste_disturbances; 2, CoV_no_smell_taste_disturbances; 3, healthCtr)",
"ID"="Number",
"salt"="How much salt do you use in your  meals?",
"bitter"="How much bitter additives do you use in your  meals?",
"sweet"="How much sugar/sweeteners do you use in your meals?",
"sour"="How much sour additives (e.g. citric acid, vinegar) do you use in your  meals?")%>%
mutate(Group=recode(Group, "1"="CoV_smell_taste_disturbances",
"2"="CoV_no_smell_taste_disturbances","3"="healthCtr"))%>%
drop_na(ID)%>%
mutate_all(~ifelse(.=="None",0, ifelse(.=="Little",1,ifelse(.=="VeryLittle", 2,
ifelse(.=="Moderately",3, ifelse(.=="Much",4, ifelse(.=="VeryMuch", 5,.)))))))%>%
mutate(salt=as.numeric(salt),
sweet=as.numeric(sweet),
bitter=as.numeric(bitter),
sour=as.numeric(sour))%>%
glimpse()
## Rows: 71
## Columns: 7
## $ Group  <chr> "CoV_smell_taste_disturbances", "CoV_smell_taste_disturbances",…
## $ ID     <chr> "GDK003", "GDK004", "GDK007", "GDK011", "GDK012", "GDK021", "GD…
## $ Sex    <chr> "M", "M", "F", "F", "F", "F", "F", "M", "F", "F", "F", "F", "M"…
## $ salt   <dbl> 3, 3, 3, 3, 3, 3, 4, 4, 2, 4, 2, 3, 4, 4, 3, 1, 3, 3, 3, 3, 4, …
## $ bitter <dbl> 3, 1, 4, 3, 3, 5, 4, 1, 1, 4, 3, 0, 3, 3, 3, 3, 2, 1, 1, 1, 4, …
## $ sweet  <dbl> 3, 5, 3, 2, 1, 1, 4, 3, 2, 2, 3, 0, 2, 2, 3, 1, 3, 0, 3, 4, 1, …
## $ sour   <dbl> 3, 3, 1, 1, 2, 0, 3, 3, 3, 3, 1, 2, 3, 3, 1, 2, 2, 2, 1, 2, 3, …
summary(df)
##     Group                ID                Sex                 salt      
##  Length:71          Length:71          Length:71          Min.   :1.000  
##  Class :character   Class :character   Class :character   1st Qu.:2.000  
##  Mode  :character   Mode  :character   Mode  :character   Median :3.000  
##                                                           Mean   :2.662  
##                                                           3rd Qu.:3.000  
##                                                           Max.   :4.000  
##      bitter          sweet            sour     
##  Min.   :0.000   Min.   :0.000   Min.   :0.00  
##  1st Qu.:1.000   1st Qu.:1.000   1st Qu.:1.00  
##  Median :3.000   Median :2.000   Median :2.00  
##  Mean   :2.507   Mean   :2.324   Mean   :2.07  
##  3rd Qu.:3.000   3rd Qu.:3.000   3rd Qu.:3.00  
##  Max.   :5.000   Max.   :5.000   Max.   :4.00
ggbetweenstats(
  data = df,
  x = Group,
  y = salt,
  title="Salt"
)

ggbetweenstats(
  data = df,
  x = Group,
  y = sweet,
  title="Sweet"
)

ggbetweenstats(
  data = df,
  x = Group,
  y = bitter,
  title="bitter"
)

ggbetweenstats(
  data = df,
  x = Group,
  y = sour,
  title="Sour"
)

# subset data
df1<-subset(df, Group !="CoV_no_smell_taste_disturbances")
ggbetweenstats(
  data = df1,
  x = Group,
  y = salt,
  title="Salt : CoV_smell_taste_disturbances vs healthCtr"
)

ggbetweenstats(
  data = df1,
  x = Group,
  y = sweet,
  title="Sweet : CoV_smell_taste_disturbances vs healthCtr"
)

ggbetweenstats(
  data = df1,
  x = Group,
  y = bitter,
  title="bitter : CoV_smell_taste_disturbances vs healthCtr"
)

ggbetweenstats(
  data = df1,
  x = Group,
  y = sour,
  title="Sour : CoV_smell_taste_disturbances vs healthCtr"
)

#salt
scheirerRayHare(salt~Group+Sex, data=df)
## 
## DV:  salt 
## Observations:  71 
## D:  0.8826626 
## MS total:  426
##           Df  Sum Sq      H p.value
## Group      2   369.9 0.9836 0.61151
## Sex        1    32.9 0.0875 0.76737
## Group:Sex  2  1930.0 5.1328 0.07681
## Residuals 65 23988.2
scheirerRayHare(salt~Group+Sex, data=subset(df, Group != "CoV_no_smell_taste_disturbances"))
## 
## DV:  salt 
## Observations:  64 
## D:  0.8804716 
## MS total:  346.6667
##           Df  Sum Sq      H p.value
## Group      1   263.4 0.8629 0.35293
## Sex        1     0.3 0.0009 0.97566
## Group:Sex  1  1381.3 4.5255 0.03339
## Residuals 60 17583.2
summary(aov(salt~Group*Sex, data=subset(df, Group != "CoV_no_smell_taste_disturbances")))
##             Df Sum Sq Mean Sq F value Pr(>F)  
## Group        1   0.67   0.673   0.733  0.395  
## Sex          1   0.00   0.003   0.003  0.956  
## Group:Sex    1   4.33   4.325   4.709  0.034 *
## Residuals   60  55.11   0.918                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#sweet
scheirerRayHare(sweet~Group+Sex, data=df)
## 
## DV:  sweet 
## Observations:  71 
## D:  0.9414655 
## MS total:  426
##           Df  Sum Sq      H p.value
## Group      2   341.5 0.8515 0.65328
## Sex        1   239.7 0.5976 0.43948
## Group:Sex  2  1500.1 3.7403 0.15410
## Residuals 65 26005.3
scheirerRayHare(sweet~Group+Sex, data=subset(df, Group != "CoV_no_smell_taste_disturbances"))
## 
## DV:  sweet 
## Observations:  64 
## D:  0.9405678 
## MS total:  346.6667
##           Df  Sum Sq      H p.value
## Group      1   260.6 0.7991 0.37136
## Sex        1    80.7 0.2474 0.61888
## Group:Sex  1  1066.2 3.2698 0.07057
## Residuals 60 19146.9
summary(aov(sweet~Group*Sex, data=subset(df, Group != "CoV_no_smell_taste_disturbances")))
##             Df Sum Sq Mean Sq F value Pr(>F)  
## Group        1   0.95   0.954   0.621  0.434  
## Sex          1   0.69   0.691   0.450  0.505  
## Group:Sex    1   5.93   5.932   3.862  0.054 .
## Residuals   60  92.17   1.536                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#bitter
scheirerRayHare(bitter~Group+Sex, data=df)
## 
## DV:  bitter 
## Observations:  71 
## D:  0.9156439 
## MS total:  426
##           Df  Sum Sq       H p.value
## Group      2   981.6 2.51663 0.28413
## Sex        1    16.5 0.04222 0.83719
## Group:Sex  2   553.4 1.41882 0.49194
## Residuals 65 25707.5
scheirerRayHare(bitter~Group+Sex, data=subset(df, Group != "CoV_no_smell_taste_disturbances"))
## 
## DV:  bitter 
## Observations:  64 
## D:  0.9123397 
## MS total:  346.6667
##           Df  Sum Sq       H p.value
## Group      1     0.3 0.00083 0.97697
## Sex        1    33.3 0.10528 0.74559
## Group:Sex  1   418.6 1.32341 0.24998
## Residuals 60 19473.6
summary(aov(bitter~Group*Sex, data=subset(df, Group != "CoV_no_smell_taste_disturbances")))
##             Df Sum Sq Mean Sq F value Pr(>F)
## Group        1   0.01  0.0135   0.009  0.926
## Sex          1   0.18  0.1831   0.120  0.731
## Group:Sex    1   1.56  1.5638   1.022  0.316
## Residuals   60  91.85  1.5308
#sour
scheirerRayHare(sour~Group+Sex, data=df)
## 
## DV:  sour 
## Observations:  71 
## D:  0.9074111 
## MS total:  426
##           Df  Sum Sq      H p.value
## Group      2   439.0 1.1356 0.56678
## Sex        1  1273.9 3.2955 0.06947
## Group:Sex  2  1031.2 2.6677 0.26346
## Residuals 65 24141.3
scheirerRayHare(sour~Group+Sex, data=subset(df, Group != "CoV_no_smell_taste_disturbances"))
## 
## DV:  sour 
## Observations:  64 
## D:  0.8957875 
## MS total:  346.6667
##           Df  Sum Sq       H p.value
## Group      1   269.3 0.86731 0.35170
## Sex        1   495.4 1.59527 0.20658
## Group:Sex  1   480.8 1.54812 0.21341
## Residuals 60 18284.3
summary(aov(sour~Group*Sex, data=subset(df, Group != "CoV_no_smell_taste_disturbances")))
##             Df Sum Sq Mean Sq F value Pr(>F)
## Group        1   1.17   1.168   1.067  0.306
## Sex          1   1.66   1.665   1.521  0.222
## Group:Sex    1   1.42   1.417   1.294  0.260
## Residuals   60  65.69   1.095

Among several advanced statistical methods and resampling techniques, Bayesian can provide more robust insights. We use Bayesian model to perform analyses of the two-way (Group * Sex) data

# salt: almost no Group difference
model<-brm(salt~Group*Sex, data=df1, family=gaussian())
# Posterior Distributions
summary(model)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: salt ~ Group * Sex 
##    Data: df1 (Number of observations: 64) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Regression Coefficients:
##                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept               2.68      0.22     2.25     3.11 1.00     3642     3130
## GrouphealthCtr          0.00      0.27    -0.55     0.55 1.00     3529     2437
## SexM                    0.82      0.53    -0.23     1.87 1.00     2251     2394
## GrouphealthCtr:SexM    -1.49      0.72    -2.90    -0.04 1.00     2621     2529
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.98      0.09     0.81     1.18 1.00     4036     2877
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# Extract posterior samples
post_samples <- posterior_samples(model)
# Plot posterior distributions of group effect
mcmc_dens(post_samples, pars = "b_GrouphealthCtr")+
  ggtitle("Salt: Posterior Distribution of Group healthCtr Effect")

# Posterior predictive checks
pp_check(model)+
  ggtitle("Salt: Posterior Predictive Check for Group Differences")

# Credible Intervals Plot
draws <- as_draws_df(model)
mcmc_trace(draws)

plot(mcmc_intervals(as_draws_df(model)))+
  ggtitle("Salt: 95% Credible Interval for Group healthCtr Effect")

# Difference in Group Means Plot 

#extracting group-level means from the posterior samples
CoV_smell_taste_disturbances_mean <- post_samples$b_Intercept
healthCtr_mean <- post_samples$b_Intercept + post_samples$b_GrouphealthCtr

#Difference in means
group_difference <- healthCtr_mean - CoV_smell_taste_disturbances_mean

#Plot the difference in means
mcmc_areas(as.data.frame(group_difference)) +
  ggtitle("Salt: Posterior Distribution of Group Differences\n (Group_healthCtr - Group_CoV_smell_taste_disturbances)")

# Interaction Plot (for Group * Sex Interaction)
new_data <- data.frame(
  Group = factor(c("CoV_smell_taste_disturbances", "healthCtr")),
  Sex = factor(c("M", "F"))
)
fitted_values <- posterior_epred(model, newdata = new_data)

# Convert to a tidy format
fitted_values_df <- add_epred_draws(new_data, model)

# Plot interaction between Group and Sex
ggplot(fitted_values_df, aes(x = Group, y = .epred, color = Sex)) +
  stat_pointinterval() +
  ggtitle("Salt: Predicted Group Differences by Sex")+theme_ggdist()

## Raincloud Plot for Group Difference
ggplot(df1, aes(x = Group, y = salt, fill = Group)) +
  stat_halfeye(adjust = .5, justification = -.2, .width = 0, point_colour = NA) +
  geom_boxplot(width = .12, outlier.color = NA, alpha = 0.5) +
  geom_jitter(width = .1, alpha = 0.3) +
  ggtitle("Salt: Raincloud Plot of Group Differences")+theme_ggdist()+
 theme(legend.position = "none")

# sweet:COVID patients reduced sugar intake
model<-brm(sweet~Group*Sex, data=df1, family=gaussian())
# Posterior Distributions
summary(model)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: sweet ~ Group * Sex 
##    Data: df1 (Number of observations: 64) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Regression Coefficients:
##                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept               1.95      0.27     1.41     2.48 1.00     3275     3120
## GrouphealthCtr          0.52      0.35    -0.15     1.20 1.00     3382     2630
## SexM                    1.28      0.67    -0.00     2.59 1.00     2399     3043
## GrouphealthCtr:SexM    -1.75      0.90    -3.48     0.03 1.00     2263     3001
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     1.26      0.11     1.06     1.51 1.00     3222     2973
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# Extract posterior samples
post_samples <- posterior_samples(model)
# Plot posterior distributions of group effect
mcmc_dens(post_samples, pars = "b_GrouphealthCtr")+
  ggtitle("Sweet: Posterior Distribution of Group healthCtr Effect")

# Posterior predictive checks
pp_check(model)+
  ggtitle("Sweet: Posterior Predictive Check for Group Differences")

# Credible Intervals Plot
draws <- as_draws_df(model)
mcmc_trace(draws)

plot(mcmc_intervals(as_draws_df(model)))+
  ggtitle("Sweet: 95% Credible Interval for Group healthCtr Effect")

# Difference in Group Means Plot 

#extracting group-level means from the posterior samples
CoV_smell_taste_disturbances_mean <- post_samples$b_Intercept
healthCtr_mean <- post_samples$b_Intercept + post_samples$b_GrouphealthCtr

#Difference in means
group_difference <- healthCtr_mean - CoV_smell_taste_disturbances_mean

#Plot the difference in means
mcmc_areas(as.data.frame(group_difference)) +
  ggtitle("Sweet: Posterior Distribution of Group Differences\n (Group_healthCtr - Group_CoV_smell_taste_disturbances)")

# Interaction Plot (for Group * Sex Interaction)
new_data <- data.frame(
  Group = factor(c("CoV_smell_taste_disturbances", "healthCtr")),
  Sex = factor(c("M", "F"))
)
fitted_values <- posterior_epred(model, newdata = new_data)

# Convert to a tidy format
fitted_values_df <- add_epred_draws(new_data, model)

# Plot interaction between Group and Sex
ggplot(fitted_values_df, aes(x = Group, y = .epred, color = Sex)) +
  stat_pointinterval() +
  ggtitle("Sweet: Predicted Group Differences by Sex")+theme_ggdist()

## Raincloud Plot for Group Difference
ggplot(df1, aes(x = Group, y = sweet, fill = Group)) +
  stat_halfeye(adjust = .5, justification = -.2, .width = 0, point_colour = NA) +
  geom_boxplot(width = .12, outlier.color = NA, alpha = 0.5) +
  geom_jitter(width = .1, alpha = 0.3) +
  ggtitle("Sweet: Raincloud Plot of Group Differences")+theme_ggdist()+
 theme(legend.position = "none")

## bitter: COVID patients increase bitter additive intake!
model<-brm(bitter~Group*Sex, data=df1, family=gaussian())
# Posterior Distributions
summary(model)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: bitter ~ Group * Sex 
##    Data: df1 (Number of observations: 64) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Regression Coefficients:
##                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept               2.67      0.27     2.14     3.20 1.00     3203     2876
## GrouphealthCtr         -0.11      0.35    -0.79     0.56 1.00     3069     2968
## SexM                   -0.66      0.70    -2.01     0.71 1.00     2411     2785
## GrouphealthCtr:SexM     0.90      0.91    -0.89     2.66 1.00     2277     2544
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     1.26      0.12     1.06     1.52 1.00     4197     2520
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# Extract posterior samples
post_samples <- posterior_samples(model)
# Plot posterior distributions of group effect
mcmc_dens(post_samples, pars = "b_GrouphealthCtr")+
  ggtitle("Bitter: Posterior Distribution of Group healthCtr Effect")

# Posterior predictive checks
pp_check(model)+
  ggtitle("Bitter: Posterior Predictive Check for Group Differences")

# Credible Intervals Plot
draws <- as_draws_df(model)
mcmc_trace(draws)

plot(mcmc_intervals(as_draws_df(model)))+
  ggtitle("Bitter: 95% Credible Interval for Group healthCtr Effect")

# Difference in Group Means Plot 

#extracting group-level means from the posterior samples
CoV_smell_taste_disturbances_mean <- post_samples$b_Intercept
healthCtr_mean <- post_samples$b_Intercept + post_samples$b_GrouphealthCtr

#Difference in means
group_difference <- healthCtr_mean - CoV_smell_taste_disturbances_mean

#Plot the difference in means
mcmc_areas(as.data.frame(group_difference)) +
  ggtitle("Bitter: Posterior Distribution of Group Differences\n (Group_healthCtr - Group_CoV_smell_taste_disturbances)")

# Interaction Plot (for Group * Sex Interaction)
new_data <- data.frame(
  Group = factor(c("CoV_smell_taste_disturbances", "healthCtr")),
  Sex = factor(c("M", "F"))
)
fitted_values <- posterior_epred(model, newdata = new_data)

# Convert to a tidy format
fitted_values_df <- add_epred_draws(new_data, model)

# Plot interaction between Group and Sex
ggplot(fitted_values_df, aes(x = Group, y = .epred, color = Sex)) +
  stat_pointinterval() +
  ggtitle("Bitter: Predicted Group Differences by Sex")+theme_ggdist()

## Raincloud Plot for Group Difference
ggplot(df1, aes(x = Group, y = bitter, fill = Group)) +
  stat_halfeye(adjust = .5, justification = -.2, .width = 0, point_colour = NA) +
  geom_boxplot(width = .12, outlier.color = NA, alpha = 0.5) +
  geom_jitter(width = .1, alpha = 0.3) +
  ggtitle("Bitter: Raincloud Plot of Group Differences")+theme_ggdist()+
 theme(legend.position = "none")

## sour COVID patients increase sour additive intake
model<-brm(sour~Group*Sex, data=df1, family=gaussian())
summary(model)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: sour ~ Group * Sex 
##    Data: df1 (Number of observations: 64) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Regression Coefficients:
##                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept               2.05      0.23     1.58     2.50 1.00     2990     2498
## GrouphealthCtr         -0.14      0.30    -0.73     0.44 1.00     2879     2438
## SexM                    0.94      0.58    -0.17     2.07 1.00     2299     2331
## GrouphealthCtr:SexM    -0.83      0.77    -2.29     0.63 1.00     2047     2605
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     1.07      0.10     0.90     1.28 1.00     4402     3014
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
post_samples <- posterior_samples(model)
# Plot posterior distributions of group effect
mcmc_dens(post_samples, pars = "b_GrouphealthCtr")+
  ggtitle("Sour: Posterior Distribution of Group healthCtr Effect")

# Posterior predictive checks
pp_check(model)+
  ggtitle("Sour: Posterior Predictive Check for Group Differences")

# Credible Intervals Plot
draws <- as_draws_df(model)
mcmc_trace(draws)

plot(mcmc_intervals(as_draws_df(model)))+
  ggtitle("Sour: 95% Credible Interval for Group healthCtr Effect")

# Difference in Group Means Plot 

#extracting group-level means from the posterior samples
CoV_smell_taste_disturbances_mean <- post_samples$b_Intercept
healthCtr_mean <- post_samples$b_Intercept + post_samples$b_GrouphealthCtr

#Difference in means
group_difference <- healthCtr_mean - CoV_smell_taste_disturbances_mean

#Plot the difference in means
mcmc_areas(as.data.frame(group_difference)) +
  ggtitle("Sour: Posterior Distribution of Group Differences\n (Group_healthCtr - Group_CoV_smell_taste_disturbances)")

# Interaction Plot (for Group * Sex Interaction)
new_data <- data.frame(
  Group = factor(c("CoV_smell_taste_disturbances", "healthCtr")),
  Sex = factor(c("M", "F"))
)
fitted_values <- posterior_epred(model, newdata = new_data)

# Convert to a tidy format
fitted_values_df <- add_epred_draws(new_data, model)

# Plot interaction between Group and Sex
ggplot(fitted_values_df, aes(x = Group, y = .epred, color = Sex)) +
  stat_pointinterval() +
  ggtitle("Sour: Predicted Group Differences by Sex")+theme_ggdist()

## Raincloud Plot for Group Difference
ggplot(df1, aes(x = Group, y = sour, fill = Group)) +
  stat_halfeye(adjust = .5, justification = -.2, .width = 0, point_colour = NA) +
  geom_boxplot(width = .12, outlier.color = NA, alpha = 0.5) +
  geom_jitter(width = .1, alpha = 0.3) +
  ggtitle("Sour: Raincloud Plot of Group Differences")+theme_ggdist()+
 theme(legend.position = "none")